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ABSTRACT 



We present the new results of the two-dimensional numerical experiments 
on the cross-sectional evolution of a twisted magnetic flux tube rising from the 
deeper solar convection zone (—20, 000 km) to the corona through the surface. 
The initial depth is ten times deeper than most of previous calculations focusing 
on the flux emergence from the uppermost convection zone. We find that the 
evolution is illustrated by the two-step process described below: the initial tube 
rises due to its buoyancy, subject to aerodynamic drag due to the external flow. 
Because of the azimuthal component of the magnetic field, the tube maintains 
its coherency and does not deform to become a vortex roll pair. When the 
flux tube approaches the photosphere and expands sufficiently, the plasma on 
the rising tube accumulates to suppress the tube's emergence. Therefore, the 
flux decelerates and extends horizontally beneath the surface. This new finding 
owes to our large scale simulation calculating simultaneously the dynamics within 
the interior as well as above the surface. As the magnetic pressure gradient 
increases around the surface, magnetic buoyancy instability is triggered locally 
and, as a result, the flux rises further into the solar corona. We also find that 
the deceleration occurs at a higher altitude than in our previous experiment 
using magnetic flux sheets (Toriumi and Yokoyama). By conducting parametric 
studies, we investigate the conditions for the two-step emergence of the rising 
flux tube: field strength > 1.5 x 10 4 G and the twist > 5.0 x 10" 4 km" 1 at 
-20, 000 km depth. 
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Introduction 



Emerging magnetic fluxes are thought to be the source of solar active regions (IParker 



1955). It is w idely accepted that the emerging flux has a form of an f2-shaped flux tube. 



Parker! ( 119751 ) studied the rise time of a flu x tube that eme rges due to magnetic buoyancy 



under the resistance of aerodynamic drag. 



Schiisslerl (Il979l ) carried out two-dimensional 



magnetohydrodynamic (MHD) simulations and found that the cross-section of the emerging 
flux tube is deformed by the surrounding flows into an "umbrella shape" with a pair of 



count er-rotati ng flux rolls, leading to fragmentation of the tube (see also 



19961 ). After 



Spruit 



Longcope et al. 



( Il98ll ) introduced the thin-flux-tube (TFT) approximation, the 



emergence of the flux tube has been studied b y using this model without considering the 



effect of the cross-sectional 



199- 



and 



Cali gari et al. 



1995 



fragmentation (e.g. 



' D'Silva fc Choudhuri 



Moreno- Insertis et al. 



Emonet fc Moreno-Insertis 



19951 ). 



1993 



Fan et al. 



Moreno-Insertis fc Emonet 



1993 



(1996) 



1998) carried out compressible MHD simulations using the 



twisted flux tube, and found that the az imuthal field can suppress the splitting of the rising 



tube into two vortex filaments (see also 



Porch 



20071 ). 



Numerical experiments on the f 
out widely in the last two decades. 



ux emergence above the solar surface have been carried 



Shibata et al. 



( 119891 ) produced a pioneering work on 



the nonlinear evolution of the undular mode of magnetic buoyancy instability (the Parker 
instability). The expansion of the magnetic fields into the corona was found to be in a 



self-similar way. iMagaral ( 1200 ll ) assumed a twisted magnetic flux tube that rises from the 
uppermost convection zone to the solar corona. When the apex of his flux tube entered the 
isothermal (i.e., convectively stable) photosphere, the tube extended horizontally around 
the surface. The rising tube did not develop a wake with a vortex roll pair, because the 
initial tube was located just beneath the surface (z = —1800 km) and the tube directly 
entered the photosphere before the wake is formed. 
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The transition of the rising flux throug h the solar surface has w i dely been discussed: th e 
thorou gh review in this matter is found in 



Moreno- Insertid (120061 ) 



Matsumoto h Shibata 



( 119921 ) carried out the first three-dimensional study on flux emergence, and, since th en, 



various three-dimensional experiments have been done (e.g . 



Archontis et al 



Galsgaard et al. 



2004 



2007 



Isobe et al. 



2005 



Toriumi et al. 



Murray et al. 



2006 



Dorch 


1999 




Fan 


2001 


Manchester et al. 


2004 



20111 ) . Recently, MHD calculations including radiative 



effect are used to reveal the interaction between flux emergence and the connective m otion 



near the surface ( ICheung et al. 



2007 



2008 



Tortosa-Andreu fc Moreno- Insertis 



20091 ). 



Toroidal flux tubes are also considered as an initia 



emergence (IHood et al. 



2009 



MacTaggart fc Hood 



cond ition to study the dynamics of flux 



2009J). 



The works introduced above can be divided into two groups: one is for the dynamics 
within the interior from ~ —200, 000km to the surface, and the other for the emergence 
from the uppermost convection zone (a few 1000 km depth) to the corona. An additional 
experiment is required on the evolution of the emerging flux tube from the deep convection 
layer to the corona through the surface for the furth er understanding of the flux emergence 



in a self-consistent manner (Abbett fc Fisher 



20031 ). In this paper, we perform a two- 



dimensional compressible MHD study on the buoyant twisted magnetic flux tube initially 
embedded deep in the convection zone [z = —20, 00 km), of which the depth is the same 



as our previous study (jToriumi fc Yokoyama 



20101 . hereafter Paper I) and is ten times 



deeper than most of previous calculations considering the emergence from the uppermost 
convection zone. 

In Paper I, we carried out numerical experiments on the emerging flux sheet, which 
decelerates in the convection zone and extends horizontally beneath the surface. When the 
magnetic flux accumulates near the surface due to the successive emergence below, the flux 
sheet is subject to the Parker instability and thus further evolution to the corona occurs 
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(the "two-step emergence" model). When the field is too weak, the rising flux cannot evolve 
further to the corona and remains within the solar interior ("failed emergence"), while, 
in the case of a strong field, the flux directly emerges to the upper atmosphere without 
showing a deceleration ("direct emergence"). 

The results of the present study are similar to the previous experiments. In the 
convection zone, the cross- s ection al ev olution of the flux tube is s i milar to the papers by 



Moreno-Insertis fc Emonet 



(1199a ) and 



Emonet Moreno-Insertid (119981 ). We found that 



as the flux tube rises, aerodynamic drag becomes more effective, since the external flow 
from the apex forms a wake be hind the main t ube. Also, the emergence above the surface 
is similar to the calculations by iMagaral (120011 ) . However, there is an important difference 



from those papers; the deceleration of the rising tube occurs when the tube is far below 
the photosphere, and the tube makes a flattened V-shaped structure widely beneath the 
surface. By conducting an analytic model which includes the effect of the mass pile-up, we 
confirm that the tube's deceleration is the consequence of the fluid accumulation between 
the rising tube and the photosphere. The magnetic flux stretched horizontally beneath the 
surface becomes magnetic buoyancy- unstable locally so that the second-step emergence to 
the corona takes place. These features resemble the previous two-step model in Paper I. We 
also carry out parameter surveys to study the dependence on the initial twist and the field 
strength. 

In the next section (§2]), we give a description of the model used in this study, while 
the results of the numerical experiments are shown in detail in Section El In Section HJ we 
present the results of the parameter study. Finally, in Section [5j we summarize and discuss 
the results. 
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Numerical Model 



2.1. Assumptions and Basic Equations 



We consider the buoyant rise of an isolated magnetic flux tube in the stratified ideal gas 
layers in the (y, z)-plane, where the z-coordinate increases upward. We solve the standard 
set of ideal MHD equations including constant gravitational acceleration g = (0,0, —go). 
The basic equations are the same as those of Paper I: 



! + v.(,v) = o, 



(1) 



d( P V) 



m +V-(pVV+pI 



BB B 

47T 87T 



pg = Q , 



(2) 



dB 

~dt 



V x (V x B) , 



(3) 



(pU + p + \ P V 2 )V + ^-E x B 

2 47T 



99 ■ V = , (4) 



and 



U = 



1 p 



7 - lp ' 



(5) 



£ = — V x B , 



(6) 



p = —pT , 

m 



(7) 



where U is the internal energy per unit mass, i" is the unit tensor, A; B is the Boltzmann 
constant, m is the mean molecular mass, and other symbols have their usual meanings. 
In this study, we carry out so-called 2.5-dimensional numerical simulation, that is, all the 
physical quantities are independent of x while the x-component of vector quantities (i.e., 



-7- 



velocity V and magnetic field B) is taken into account. We assume a ratio of specific heats, 
7 = 5/3. 

The normalizing units of length, velocity, time, and density in the simulations are 
H , C s o, t = Hq/C s o, and po, respectively, where H = k^To/lmgo) is the pressure 
scale height, C s0 the sound speed, and p the density at the photosphere. The gas 
pressure, temperature, and magnetic field strength are normalized by the combinations 
of the units above, i.e., po = poC 2 , T = rnC 2 / (7&b), and B = (poC^) 1 ^ 2 , respectively. 
The gravity is given as go = C 2 / (jHo) by definition. For the comparison of numerical 
results with observations, we use H Q = 200 km, C s0 = 8 kms" 1 , r = H /C s0 = 25 s, 
and po = 1-4 x 10~ 7 g cm -3 , which are typical values for the solar photosphere and 
chromosphere. Then, p = 9.0 x 10 4 dyn cm -2 , T = 4000 K, and B = 300 G. 



2.2. Initial Conditions 



The initial background stratification consists of three regions (see e.g. iNozawa et al. 



19921 ): an adiabatically stratified convection zone, a cool isothermal photosphere/chromosphere, 
and a hot isothermal solar corona. We take z = to be the base height of the photosphere, 
and the initial temperature distribution of the photosphere/chromosphere and the corona 
(z > 0) is assumed to be 

T = T s {z) = T ph + (T cor - T ph ){tanh [{z - z cm )/w tr ] + l}/2 , (8) 

where T p h = T and T cor = 100T are the respective temperatures in the photo- 
sphere/chromosphere and the corona, z COT = 10H is the base of the corona, and w tT = 0.5H 
is the temperature scale height of the transition region. The initial temperature distribution 
in the convection zone (z < 0) is described as 

dT 



T = T s (z) = T ph - z 



dz 



(9) 

ad 
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where 



dT 

dz 



7-1 mg Q 



(10) 



ad 7 k B 

is the adiabatic temperature gradient, i.e., the initial temperature distribution in the 
convection zone is adiabatic. On the basis of the temperature distribution above, the initial 
pressure and density profiles are defined by the equation of static pressure balance: 

dp s (z) 



dz 



+ Ps(z)g = 0. 



The initial magnetic flux tube is embedded in the solar interior at z — — lOOf^o 
-20, 000 km. The longi tudinal an d azimuthal components of the flux tube are 



described as follows (see 



Fan 



2001 



Archontis et al. 



2004 : 



Murray et al 



200J): for 



r=[(y~ 2/tube) 2 + (z- -2 tu bc) 2 ] 1/2 > 

B x (r) = B tnhe exp 
and 

B<j,{r) = qrB x (r) 



p2 

-"-tube 



(12) 



(13) 



where (y t ube, ^tube) = (0, — 100if ) is the tube center, -R tu be is the tube radius, q is the twist 
parameter denoting the angular rate of field lines rotating around the tube's axis per unit 
length, and -Btube is the field strength at the center. The horizontal and vertical components 
of the azimuthal field are defined as 



Z ^tube 



(14) 



and 



B z (y, z) = B^r) 



y - z/tubc 



(15) 
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respectively. We use -Rtube = 5H = 1000 km throughout the paper. For the typical case 
(case 1), we take -Btube = 50B = 1.5 x 10 4 G, so the total longitudinal magnetic flux is 
$ x = 4.7 x 10 20 Mx. The twist parameter is q = 0.1/H for case 1. The pressure inside the 
tube is defined as 



for the pressure balance. The temperature is kept unchanged, i.e., thermal balance is 
sustained, T\(z) = T s (z). Since the density inside the tube is smaller than that outside, 
the flux tube will buoyantly rise through the convection zone. The initial temperature, 
density, and pressure profiles of the background stratification, and the total field strength 



B =[B 2 X + BV\ l/2 = [B 2 X + B 2 + B 2 } 1 / 2 of case 1 along y = are shown in Figure [U 



The domain of the simulation box is (y m i n < y < y max ) and (z min < z < z max ), 
where y min = -400H , y max = 400H , z min = -150H , and z max = 250H . The total 
size of the simulation box is 160 Mm x 80 Mm. Periodic boundaries are assumed 
for y = y m i n and y = 2/ max , symmetric boundaries for z = 2 m ; n and z = z max . A 
wave-damping region is attached near the top boundary. The total number of grid points 
is (N y x N z ) = (1536 x 1920), and the mesh sizes are Ay = 0.52H and Az = 0.21if Q , 
both of which are uniform. We use the modified Lax-Wendroff scheme version of the 
CANS (Coordinated Astronomical Numerical Software) code (see Paper I). To study the 
dependence on the twist and the field strength of the flux tube, we change the values for 



Pi(y, z) = Ps(z) + 6p cxc (r) , 



(16) 



where Sp exc (< 0) is the pressure excess described as 




(17) 



2.3. Boundary Conditions and Numerical Procedures 
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the twist parameter q or the axial field strength B tu ^ e . The cases we use in this study are 
summarized in Table [TJ 



Results 



Here, we show the results of the typical model (case 1) in which the initial twisted 
flux tube with q = 0.1/ Ho rises through the convection zone to the corona. Th e evolution 



within the interior resemb l es the studies by 



Moreno-Insertis & Emonet 



(119961 ) and 



Emonet fc Moreno- Insert id (119981 ). and the evolution to the upper atmosphere is similar to 
Magaral (|200ll ). However, the present result is not a simple sum of them; the rising tube 
decelerates within the interior far below the surface to build a widely flattened magnetic 
structure. 

Figure [2] displays the evolution of the density profile (color contour) in a vertical plane 
normal to the tube axis, while the solid lines and arrows indicate the magnetic field lines 
and the velocity vectors, respectively. Figure [2]^a) presents the initial state; the flux tube 
is located at z/H = —100 with a circular shape, and begins to emerge because of the 



buoyancy due to its smaller density relative to the surroundings. Figure 3(a) shows the 



height-time relations at the top (solid line), the center (dotted line), and the bottom (dashed 



line) of the emerging tube, while the velocities at these points are presented in Figure 3(b) 



In the initial phase (0 < t/r < 100), the flux tube rises due to magnetic buoyancy and 
the tube's cross-section almost keeps its original circular shape. Thus, the tube moves up 
with a constant acceleration rate. In the next phase (100 < t/ro < 500), the aerodynamic 
drag grows to counteract buoyancy, since the external flow around the tube's cross-section 
develops a wake behind the main tube (Figure [2](c)), and the acceleration of the rise velocity 
is reduced. The azimuthal field of the flux tube resists this deformation and prevents the 
tube from fragmenting into a counter-rotating vortex pair (the so-called "umbrella shape" 
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111 



Schiisslerl (119791 ) ). During the emergence in the convection zone, the rise velocity of the 



tube center fluctuates due to the internal torsional oscillation (Figure 3(b) see §3. 1 1) . As the 
flux tube reaches close to the photosphere, it expands horizontally (Figure Efe)) and its rise 
speed decreases after t/r = 500. Eventually, the subsurface flux extends widely in the range 
of -100 < y/H < 100 at t/r = 800 (Figure EJ^f)). When the top of the subphotospheric 
flux becomes unstable to the magnetic buoyancy instability, further evolution to the corona 
breaks out (Figures [2](f )-(i) ; see §3.3p . The consequent coronal loop at the time t/r = 980 
has its width ~ A00H = 80, 000 km and height ~ 200# = 40, 000 km. The picture of this 
evolution resembles our previous "two-step emergence" model (Paper I). In the following 
sections, we will describe some more detailed results, especially on the magnetic and the 
flow field within the convection zone, the comparison with the analytic model, and the 
second-step evolution starting at the photosphere. 



3.1. Magnetic Fields and Vorticity 



To study the flux tube's emergence in the solar interior, we show the variation of the 
field configuration with time in Figure HJ The color contour indicates the distribution of the 
longitudinal field {B x ) while the strength of the azimuthal field {B^) is overplotted with 
solid lines. In Figure |5l we plot the x-component of the vorticity vector (uj x ), where 



dV z dV„ 



dy dz 

and white and black colors indicate positive (counter-clockwise) and negative vorticity 
(clockwise), respectively. These figures (Figures B and [5]) look similar to those o : 



previous experiments within the convection zone (iMoreno-Insertis fc Emonet 



1996 



Emonet &: Moreno- Insertis 



19981 ). Figures H](a) and[5](a) present the initial state: the flux 



tube has a circular shape and there is no vorticity within the domain. As the tube begins 
to emerge, the field strength reduces and the tube expands, since the background density 
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decreases with height (Figures H(b) -(d)). Also, the external flow from the apex to the flanks 
of the tube creates negative and positive vortex sheets in the right and left side of the 
main tube, respectively (Figures E^b)-(d)). Therefore, the tube is deformed and the wake 
behind the main tube grows with a pair of counter-rotating vortex rolls. Because of the 
azimuthal fields, however, the tube is not entirely fragmented into an "umbrella shape." 
As a consequence, the rising tube feels aerodynamic drag, and, thus, the rise velocity 
levels off at this time. The longitudinal field component keeps its strength in the core 
of the vortex rolls as well as in the main tube (see Figure Hie)). Behind the main tube, 
in —100 < z/Hq < —60 in Figure IH(f), there is a weakly magnetized tail with non-zero 
vorticity (Figure E^f)). 

During this phase, the rise speed of the tube center reveals oscillation as seen in Figure 



3(b) This is a result of torsion al oscillation due to the differentia l buoyancy caused by the 



magnetic field distribution (e.g. 



Moreno-Insertis fc Emonet 



19961 ). Initially, the flux tube 



has its maximum field strength at the center and thus the central region rises faster than the 
periphery. Therefore, the azimuthal field around the apex is compressed and strengthened 
to fortify the magnetic tension force, which inhibits the tube from a distortion into vortex 
rolls. The tension force decelerates the center of the tube and, as a result, induces the 
internal torsional oscillation. The period of this oscillation is observed to be consistent with 
the time of the azimuthal Alfven speed traveling across the tube's diameter. 

At around t/r = 500, the rise motion begins to decelerate and the front of the tube 
expands horizontally, since the tube is close enough to the surface. The plasma above 
the rising tube cannot move through the isothermal (i.e., convectively stable) photosphere 
so that the fluid is compressed and piles up between the front of the tube and the 
photosphere. The fluid between them suppresses the rising motion of the flux below, and 
thus, the flux extends sideways. The horizontal expansion caused by the photosphere was 
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previously reported by iMagaral (120011 ). His flux tube was not so deformed as to develop 
a wake including vortex rolls and directly entered the photosphere mostly keeping its 
original cylindrical shape, because the initial tube was located just beneath the surface 
(ztube = —1800 km) with a strong field strength (-B tu be = 7760 G). In our case, the initial 
tube was embedded so deep in the convection zone (z tu be = —20, 000 km) that the wake 
behind the main tube was formed and the drag force became more effective. The velocity 
shear between the horizontal flow above the tube and the photosphere forms vortex sheets 
at around z/Hq ~ (Figure [5(e); the positive sheet in the right half, and the negative in 
the left). Finally, at t/r = 700, the flux tube has a V-shaped structure and the magnetic 
field extends widely beneath the surface (Figure H](f)). The vortex sheet at the tube front 
collides against the subphotospheric sheet (Figure [5(f)). In this stage, the azimuthal field 



beneath the surface is muc h stronger t 



expansion of the flux tube ([Parker 



r an th e longitudinal field because of the radial 



1974 



19791). 



3.2. Comparison with an Analytic Model 



Fan et al.l (119981 ) compared their numerical results of the emerging twisted flux tube 
obtained by 2.5D anelastic M HD experiments with those of the thin-flux-tube model. 



According to 



Fan et al. 



(jl998l ). the motion of the cylindrical flux tube rising by its magnetic 



buoyancy is described as 



Ip- 



dt 



-Apg - Ci 



p\V z \V z 

' 7TR 



(19) 



where I is the enhanced inertia factor (~ 2), C*d is the drag co-efficient of order unity, and 
Ap = pi — p is the density difference between the flux tube (p ; ) and the external medium 
(p). From Equation ( !T9|) . one can calculate the rising velocity of the tube V z {t) and the 
height of the tube center z{t). We consider the effect of the mass pile-up between the flux 
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tube and the photosphere and replace Ap with an additional effect as 

Ap' = Ap + FAp acm 

= Pi-p + F[p acm (0) -p acm (t)] . (20) 

Here, 

1 r z Ph 

Pacm(0) = -r Ps (Q d( (21) 

is the plasma on the flux tube in the initial state (z ph = 0, z(0) = z t ube = —100H Q , and 
p s (z) is given by Equation (TTTjl ). and 

Pacm(t) = 77, / Ps(C) <*C (22) 

^ P h — ^ItJ J z (t) 

is the background plasma at a given time. That is, Ap acm = p acm (0) — p acm (£) corresponds to 
a density that would have been accumulated by an ideal rising sheet extending horizontally. 
To consider the effect of the tube's shape and the draining of the plasma from the apex, we 
multiply it by a factor F(< 1), which depends on the field geometry and its initial depth 
from the surface (-Rtube and z tu be )• In the present study, we vary F to fit the analytic 
model to the obtained data and assume it as a constant. The relation between F and -Rtube 
and -Ztube requires much work, which we shall leave for future research. Applying Ap' and 
dividing Equation (Tl9|) by p, we get 

at p ttH 

If pressure balance, mass and flux conservation, and adiabatic evolution are assumed, the 
buoyancy and the radius of the model tube can be defined as 



Ap' / Ap 
-9 = 9 



P VP/ 



(T + l)-z/H ir - 

(r + 1) - z tuhe /H 



+ F^E^L , (24) 
P 



R(z) = R{z = z t 



ubc J 



(T+l)-z/H 

(r + 1) - z tuhe /H 



r/2 



(25) 
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where r = 1/(7 — 1). Here we use their initial values obtained as 



Ap 



^ ^tubc 



and 



R(Z — z t-abe) 



sec 



Fan et al. 



ff(-&p/p)gB x dydz 
ff B x dydz 

SJ[(y - Z/tubc) 2 + {z- z tuhc ) 2 } B x dydz 1 1/J 
ff B x dydz 



(26) 



(27) 



19981 1 . We can calculate the time evolution of the model tube by integrating 



Equation (T53 



In Figure 6(a), the variation of the rise speed at the tube center i s presented in 



comparison with the analytic model (the dotted line is for the model by 



Fan et al 



(1998) 



using Ap and the solid line for our model using Ap'), while the heights of the tube centers 



of the numerical and the analytic model are indicated in Figure 6(b) We use Cd = 2.0 



and F = 0.001. In Figure 6(a), the solid curve levels off after £/tq = 100. T his means that 



the wake develops from this time and the drag force becomes more efficient. iParkerl ( 119751 ) 
analytically calculated the terminal velocity of the rising. Considering Equation ( Fl9l = 
and using Equations ( 1261) and ( 1271) . the terminal velocity can be obtained as 



V, 



Ap ttR 

9-7T- 

p C D 



ll/2 



0.14 C s 



(2? 



which is rather a good estimation. In the earlier phase {t/r < 400), the torsional oscillation 



due to the differential buoyancy is seen in Figure 6(a), which was mentioned in §3.1 



After t/r = 500, the rising flux tube decelerates as the mass on the tube is compressed 
and piles up, since the mass cannot persist through the convectively stable surface. This 
deceleration indicates that, before £/to = 500, the cross-sectional evolution as a whole can 
be regarded as radial, although the tube suffers aerodynamic deformation; after that time, 
however, the mass pile-up becomes much more efficient so that the tube decelerates and the 
apex expands horizontally to become V-shaped (see the solid and the dotted line in Figure 
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Figure [7(a) shows the density accumulation [p(t) — p s }/ Po, where p s is the initial 
background density profile (see Equation ffTT]) ). and the field lines with velocity vectors 
at the time t/ro = 600. It reveals that, because of the relative mass draining from the 
apex of the tube to the flanks, the density piles up in front of the rising tube around 
(y/Ho, z/Hq) = (0, —12) to form a boundary layer with a finite width (~ 40-f^o)- in the 
imaginary sheet case, however, the fluid would piles up on the rising sheet as a boundary 
layer with an infinite width, since there is no draining (see Paper I). Therefore, the factor 
F (the ratio of the actual accumulation to the imaginary one) becomes relatively small. 

To show that the contribution of ram pressure to the pressure excess (which is related 
to density excess) at the tube front is small enough, in Figure [TJb) we plot vertical 
distributions of ram pressure p(t)AV z (t)/p , where AV z (t) = V z (t) — V z {t,z = —1AH ) is 
the flow velocity relative to the rising tube, pressure excess [pit) — p s ]/po, and the total 
field strength \B\/B Q at t/r = 600 along the symmetric axis y/H = 0. In front of the 
rising tube {z/Hq = —12), pressure excess reveals a hump AP/po, which is indicated by 
an arrow in the middle of Figure UJb). The corresponding ram pressure p(t)AV z (t)/po is 
also indicated with arrows in the bottom of the figure. From this figure, we can see that 
the effect of the ram pressure p(t)AV z (t)/po ~ 0.002 is much smaller than the hump of the 
pressure excess AP/p ~ 1. That is, the contribution of the ram pressure to the density 
accumulation in front of the tube is negligible. 

From Figure [7Ja), we can estimate the order of F. The actual mass pile- up along the 
symmetric axis y/H Q = is 



where ^ ap cx(600r ) = — 12if is the height of the apex of the tube at t = 600r (see Figure 
UJ). On the other hand, the density enhancement in front of the imaginary emerging sheet 




0.73p , 



(29) 



-17- 



at this time can be estimated from the background profile p s {z) as 



Pacm (0) Pacm 

(600r 

1 f Z ph 



96.1/no, 



Ps(C) d( 



z p h - ^(600r ) 7 2 ( 6 oot ) 



Ps(C) 



(30) 



where 2:(600r ) = —25H is the height of the tube's center at t = 600r (see Figure 6(b)). 
Therefore, we can see that the order of the factor F is 



F = = = 0.007 = O(10- 2 ) 



(31) 



Ap acm 96.1p 

Here, we consider the density accumulation only along the axis y/H = in Equation (|29|) . 
The factor F would be of order 10~ 3 if we take it into account that the accumulation in the 
neighboring region is less than that along y/Ho = due to mass draining. 

By varying F to fit the analytic model to the obtained data and plotting the density 
pile-up in the interior, we can conclude that th e decelerat i on of the t ube is caused by the 



Fan et al. 



(119981 ) and 



Cheung et al. 



(120061) 



accumulation of the plasma ahead of the tube, 
also reported this deceleration. However, their calculations did not include the convectively 
stable photosphere but assumed the non-penetrating or closed top boundaries, which are 
reasonable approximations for the photosphere compared to our simulation. 



3.3. Further Evolution to the Corona 

After t/ro = 800, the second-step emergence from the surface to the corona occurs 



due to the magnetic buoyancy instability (Figures [2(f)- (i)). Figure 8(a) shows the vertical 
distribution of magnetic and gas pressure, and density along the symmetric axis y/H = 
at t/ro = 800, i.e., just before the secondary evolution starts. As can be seen from this 
figure, the site of the second-step evolution has a top heavy structure, and there is a relative 
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"pres sure hill," which is consistent with preceding studies ( IMagaral l200ll ; lArchontis et al. 
20041) . 



To confirm that further evolution is caused by t he ma gnetic buoyancy instability, we 



conduct the same analysis as in Paper I. 



Newcombl (119611 ) revealed the criterion for the 



magnetic buoyancy instability of a flat magnetized atmosphere is 

_ dp p 2 g 

dz 7p 

From this relation, we define the index 

, dp p 2 g 

i> = , 

oz 7p 



(32) 



(33) 



where the area with negative ip is magnetic buoyancy-unstable. We show the ip distribution 



in Figure 8(b) As can be seen from Figure 8(b), the site of the second-step emergence 
(y/Ho = z/H = 0) is ip < 0, i.e., the criterion ( 132]) for the instability is satisfied. Thus, 
we can conclude that the further evolution to the corona occurs because of the magnetic 
buoyancy instability. At this site, the magnetic field intensity amounts to ~1 kG, and the 
plasma beta is (3 ~ 2. 

In the upper atmosphere above the surface, the azimuthal component of the magnetic 
field is very dominant. The sec ond-step nonlinear ev olution develops in a self-similar way; 



the expansion law was given by 



Shibata et al. 



(119891 ). It can be described as follows: 



V z /C s0 = az/H 



(34) 



p oc z , 



(35) 



and 



(36) 
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where a ~ 0.05 when plasma-/? is 0.5 — 2.0. The above relations are plotted in Figure|9] for 



t/r = 850, 870, and 890. In Figure 9(a), we use a = 0.04, which is appropriate because 



plasma-/? was ~ 2 at the point of the further evolution (z ~ at the time t/r = 800). 

In the final phase, at £/ro = 980, the coronal loop is formed with 400ifo — 80, 000 km in 
width and 200H$ — 40, 000 km in height (see Figure E(i)). The size of the loop is similar to 
that of the case with a flux sheet (Paper I). The photospheric field strength is 400 — 700 G 
at this time. 



4. Dependence on the Initial Twist and the Field Strength 



We carry out five additional calculations to investigate the dependence of the 
tube's evolution on the twist and the field strength. These runs are summarized 
in Table [D Parametric su r veys o n the twist and the field strength have done by 



Moreno-Insertis fc Emonet 



Murray et al. 



(2006 



1996) 



Emonet fc Moreno-Insertis ( 



Murray k Hood! (120081 ). and 



Toriumi et al. 



998) 



Magaral fl2001[ ). 



(120111) • They found that 



the tube rises faster as the twist and the field becomes stronger. 

Figure [TU] shows the result of the parameter survey on the tube twist. In this figure, 
height-time relations of the tube tops are shown. In cases of a weaker twist (q = 0.01/ H 
and 0.05/ H ), the flux tube shows "failed emergence." That is, the rising flux cannot pass 
through the photosphere. The criterion for this failure is q < 0.05/ Hq. In Paper I, we found 
that, when the initial field strength or the total magnetic flux is too weak, the resulting 
development also reveals "failed" evolution. Figure [TT] shows the flux tubes with different 
initial twist when each tube center reaches z/Hq = —50. From left to right, the twist 
parameter decreases. The upper panels present the flow field relative to the tube apex, 
V re \ = V — V apex , and the corresponding equipartition line. The equipartition line indicates 
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the tube boundary where the kinetic energy density of the relative flow e^n = pKei/2 
equals the energy density of the azimuthal field e mag = B?/(8tt). The lower panels show the 
longitudinal field strength B x and the flow field V. As can be seen from Figure [TT], with a 
decreasing twist, the main tube is deformed by the external flow and the wake develops. In 
cases of a weaker twist, the counter-rotating vortices contain a large portion of the magnetic 
flux. Thus, when the tube reaches the surface, the apex of the tube cannot hold sufficient 
flux for the further evolution, which leads to "failed emergence." The situation can also 
be explained as follows; when the initial twist is weak, the flux tube expands and extends 
very widely beneath the surface. Thus, the magnetic buoyancy falls short of compressing 
the flux to satisfy the criterion (|32|) for the second-step emergence. In Figure [121 we show 
the distribution near the surface for the case q = 0.05/ H at the time t/r = 1000. Figure 



12(a) plots the magnetic pressure, the gas pressure, and the density along the symmetric 
axis y/H = 0. The subsurface magnetic field is not fortified enough so that the magnetic 
pressure is much smaller than the plasma pressure. Therefore, the index ip (see §3.31) is not 



negative in almost the whole area shown in Figure 12(b) which means that the further 
evolution to the solar corona due to the magnetic buoyancy instability does not take place 
in this weaker twist case. 

Figure [13] shows the height-time relations for the cases with -B tu bc = 67£>o, 50B , 
and 33i?o, which are indicated with dashed, solid, and dash-dotted lines. We can 
conclude from this figure that the rise speed is faster with the stronger field. As is 
the case with the dependence on the twist, we can say that, if the initial field is too 
weak, the rise speed is slower and the flux tube cannot pass through the surface to rise 
further into the corona. The general tendency of the parametric studies varying the 



experiments f 


More 


10-Insertis & Emonet 


1996 


; 


5monet & Moreno-I 


nsertis 


1998: 


Mascara 


2001: 


Murrav et al. 


2006; 


Murrav & Hoo J 


2008 


Toriumi et al. 


2011 
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5. Summary and Discussion 

The numerical experiments studied in this paper are on the cross-sectional evolution 
of the emerging twisted flux tube. In this section, we summarize the calculations presented 
above and discuss the results. The discussions are mainly in connection with Paper I. We 
can predict the conditions for three-dimensional experiments by considering the results 
of the present and previous two-dimensional studies. So, in this section, we discuss the 
two-dimensional results in connection with our future three-dimensional experiments. 

In this paper, we have studied the dynamical evolution of the twisted flux tube that 
emerged from the deep convection zone. For the typical case, the initial flux tube is located 
at z = -100# = -20, 000 km with the axial field strength B tuhc = 50B = 1.5 x 10 4 G, the 
tube radius -Rtube = 5Hq = 1000 km, and the total magnetic flux $ = 4.7 x 10 20 Mx. The 
tube is initially buoyant and, thus, begins to rise through the convection zone. Halfway to 
the surface, after t = 100r = 2500 s, the periphery of the tube is peeled away to develop 
a wake, composed of vortex rolls and a long-drawn tail, due to the external flow from the 
apex to the flanks. As a result, the aerodynamic drag becomes more effective. However, the 
azimuthal field prevents the flux tube from being fragmented into a pair of counter-rotating 
vortices moving away from each other, so that the expansion as a whole can be said to 
be radial. The picture of the emergence with in the convective layer is similar to previous 



papers (e.g. 



Moreno-Insertis fc Emonet 



19961 ). 



The essential difference from these preceding studies is the effect of the convectively 
stable photosphere; the apex of the tube expanding horizontally and the rising motion 
turning into the deceleration after t = 500ro = 1.3 x 10 4 s, because the rising tube comes 
close to the photosphere. The deceleration occurs at ten times deeper than previously 



reported by 



Magaral (120011 ); the deceleration depth is ~ —500 km in Niagara's case, while it 



is ~ —5000 km in our case, which shows that the flux tube slows down before the tube itself 
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reaches the surface. There seems to be another mechanism for the slowdown other than 
the tube itself entering the photosphere. The plasma on the flux tube piles up between the 
flux and the convectively stable surface, and depresses the flux below. We confirmed this 
effect by comparing the numerical result with an extended analytic model (see Section I3T2]) . 
As a result, the V-shaped emerging flux extends widely beneath the surface (~ 40, 000 km 
in width). In this stage, the azimuthal field component is dominant, since the twist of 
the flux tube increases as the tube rises. When the horizontal flux at the photosphere 
becomes unstable to the magnetic buoyancy instability at t — 800to = 2.0 x 10 4 s, the 
second-step emergence to the corona t akes place ("tw o-step emergence"). The evolution 



above the surface is similar to that of 



described by the expansion law given by 



Maearal (1200 lh. 



Shibata et al. 



'he nonlinear evolution is well 



( 119891 ) . The consequent coronal 



loop at t = 980r = 2.5 x 10 4 s has its width ~ A00H = 80, 000 km and height 
~ 200#o = 40, 000 km, and the photospheric field is (1.3 - 2.3) J B = 400 - 700 G. 

By performing parameter studies, we found that the initial twist necessary for the further 
evolution is q > 0.05/ H = 2.5 x 10~ 4 km -1 when the axial field -B tu be = 50-B = 1.5 x 10 4 
G. If q is less than this value, the emerging tube remains within the convection zone 
("failed emergence"). We also found that the flux with a weaker field strength shows the 
failed emergence. When the tube twist is q = 0.1/ Hq = 5.0 x 10" 4 km" 1 , the flux with 
-Btube = 33i?o = 1.0 x 10 4 G cannot rise further above the surface. 

In Paper I, we carried out numerical experiments on the two-dimensional undular 
evolution of the flux sheet from the deep convection zone (z = —20,000 km). We found 
that the emerging flux sheet also decelerates in the convection zone. However, deceleration 
depth is deeper than that of the present study. Table [2] summarizes the characteristic values 
of both types of emergence. As seen from Table El the initial conditions (the field strength, 
the total magnetic flux, the sheet thickness and the tube radius, and the depth) are similar 
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to each other, while the obtained values concerning the evolutions in the convection zone 
(the arrival time at the photosphere and the deceleration depth) are different. That is, the 
twisted flux tube rises faster than the flux sheet and the deceleration occurs at a higher 
altitude. The difference between the two cases comes from the geometry of the emerging 
fluxes, i.e., a sheet or a tube. The plasma on the emerging flux tube can flow around the 
tube's cross-section so that the tube rises faster. The fluid on the flux sheet, however, 
drains only along the magnetic field lines to both troughs, so that the plasma suppresses 
and slows down the rising sheet below halfway to the surface. As for the twisted flux tube, 
the deceleration occurs in the later phase when the tube itself approaches the surface and 
begins to expand horizontally (see §3.11 and §3.2p . It should be noted that, in the initial 
state, the flux sheet is in hydrostatic equilibrium while the flux tube is in mechanically 
non-equilibrium. In addition, the mass draining along the longitudinal field line is not 
considered in the present tube calculation. Three-dimensional experiments are required to 
investigate these effects. 

By conducting the parameter survey, we found that the emerging flux tube cannot rise 
further if the initial twist is too weak ("failed" case), because the weak azimuthal field 
cannot hold the tube's coherency and thus cannot hold the field intensity necessary for the 
second-step emergence. Such a flux remains below the surface and floats around the solar 
interior. The condition for the "two-step emerging" flux tube that yields a realistic active 
region at the surface is B tuhc > 50B = 1.5 x 10 4 G with q > 0.1/ H = 5.0 x 10" 4 km -1 
at z = — 100-fiTo — —20,000 km. In Paper I, it was also found that there is a threshold 
of the flux sheet strength and the total flux for further emergence. That is to say, the 
emerging flux with a weak twist also shows "failed emergence" as well as the flux with a 
weak field strength and an insufficient total flux. If we put these results obtained from 
the two experiments together, the preferable initial conditions for the three-dimensional 
calculation using the twisted flux tube from the same depth (z = —20, 000 km) are the field 



strength ~ a few 10 4 G, the total magnetic flux of ~ 10 21 — 10 22 Mx, and sufficient twist 
(q > 0.05///o = 2.5 x 1CU 4 km -1 ). In future three-dimensional simulations, we will take 
into account the above mentioned conditions. 

The comparison with the results by thin-flux-tube experiments (see JQ) has been 
discussed in Paper I. In Figure 10 of Paper I, the typical case of the present calculation 
(Case 1 with 1.5 x 1 4 G and 4.7 x 10 20 Mx) loc ates in the middle of the area for the 



two-step emergence. 



Moreno-Insertis et al 



(119951 ) found that the magnetic flux tube with 



weaker field "explodes" within the interior and never reaches the surface; when the pressure 
gap between inside and outside the initial flux tube is too small, i.e., when the initial 
magnetic field is too weak, the tube will be collapsed at a certain height, since the pressure 
gap decreases. In Paper I, we found that the cases showing two-step emergence would have 
survived the "explosion" during their ascents through the interior. Therefore, we can see 
that the typical tube (the present Case 1) also would have emerged through the convection 
zone without suffering explosion. 



Numerical computations were carried out on NEC SX-9 at the Center for Computational 
Astrophysics, CfCA, of the National Astronomical Observatory of Japan, and on M System 
(Fujitsu FX1) and V System (NEC SX-9) of JAXA Supercomputer System. The page 
charge of this paper is partly supported by CfCA. The authors thank Drs. H. Isobe and K. 
Shibata of Kyoto University and Dr. Y. Fan of the High Altitude Observatory, the National 
Center for Atmospheric Research. We are grateful to Mr. David Ward, the GCOE program 
instructor of the University of Tokyo for proofreading/editing assistance. 
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Fig. I. — One-dimensional (z-) distributions of the initial background density (solid line), 
pressure (dotted line), and temperature (dashed line). The total magnetic field strength 
B = [B% + B|] 1 I 2 of case 1 along the vertical axis y = is overplotted with a dashed-dotted 
line. 
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Fig. 2. — Time-evolution of the "two-step emergence" of the twisted flux tube (typical 
model), (a) t/r = 0; (b) t/r = 200; (c) t/r = 400; (d) t/r = 550; (e) t/r = 700; (f) 
t/ To = 800; (g) t/r = 870; (h) t/r = 900; (i) t/r = 980. Logarithmic density profiles 
(log 10 (p/po)) are indicated by color contours, while magnetic field lines and velocity vectors 
are overplotted with black lines and arrows, respectively. This figure is also available as an 
avi animation in the electronic edition of the Astrophysical Journal. 
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Fig. 3. — (a) Height-time relations at the top (solid line), the center (dotted line), and the 
bottom (dashed line) of the flux tube, (b) : Gas velocities at these three points. 



31 




-60 -40 -20 20 40 60 -60 -40 -20 20 40 60 




10 20 30 40 50 



Fig. 4. — Time-evolution of the magnetic fields from t/ro = to £/to = 700. Color con- 
tour displays the longitudinal component B x /Bq, while the azimuthal component B^/Bq is 
overplotted with solid lines. 
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Fig. 5. — Same as Figure Hlbut for the x-component of the vorticity. White indicates positive 
(counter-clockwise) vorticity, while black is negative (clockwise). 
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Fig. 6.- 



The variation of rise velocity of the rising flux tube with time. Diamonds 
indicate the result of the numerical si mulation. Solid line shows our analytic model, while 
dotted line represents the model by 



Fanetal 



(1l998l ). (b): Height-time relation of the 



numerical results (diamonds) and the analytic model (solid line). We use Cd = 2.0 and 
F = 0.001. 
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Fig. 7. — (a): The density accumulation [p(t) — p s ]/po, where p s is the background den- 
sity profile, and the field lines (contours) with velocity vectors (white arrows) at the time 
t/r = 600 are shown. As the tube rises, the mass piles up in front of the flux tube around 
(y/H ,z/H ) = (0,-12). (b): Vertical profiles of the ram pressure p{t)AV^(t)/p (solid 
line), the pressure excess [p(t) — p s ]/po (dotted line), and the total field strength \B\/B 
(dashed line) along the symmetric axis y/H = 0. The hump of the pressure excess AP/p 
is indicated by an arrow in the middle of the figure, while the corresponding ram pressure is 
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WH 

Fig. 8. — (a): Vertical distributions of the magnetic pressure (solid line), gas pressure 
(dotted line), and gas density (dashed line) along the axis y/H = at t/r = 800. (b): 
Two-dimensional map of the index ip = —dp/dz — p 2 go/{lp)- Color contour indicates ip, 
while azimuthal magnetic field lines are overplotted with solid lines. The area where index 
ip < is subject to the magnetic buoyancy instability. 
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Fig. 9. — (a): Distribution of the upward velocity along the vertical axis y/H = 0. Dotted 
dashed, and dash-dotted lines indicate the distribution at t/r Q = 850, 870, an d 890, respec- 



tively The solid line shows the theoretical velocity-height relation according to 



Shibata et al. 



( I1989I ). (b): Distribution of the gas density, (c): Distribution of the horizontal component 



of the magnetic field. 
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Fig. 10. — Height-time relation of the top of the flux tube. Dashed, solid, dash-dotted, and 
long dashed lines represent cases for q = 0.2/ H , 0.1 /H , 0.05 /H , and 0.01/ H , respectively. 
Dotted line indicates the photospheric level. 
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Fig. 11. — Cross-section of the flux tubes with different initial twist when their centers are 
at z/Hq = —50. From left to right, the twist parameters and their corresponding times are 
q = 0.2/ H , 0.1/H Q , 0.05/H , and 0M/H , and t/r = 460, 410, 390, and 450, respectively. 
Top: The velocity relative to the apex of the tube and the equipartition line (see text for 
details). Bottom: The longitudinal magnetic field strength and the flow fields. 
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Fig. 12. — The plots for the case with q = 0.05/H at t/r = 1000. (a): Vertical distributions 
of the magnetic pressure (solid line), gas pressure (dotted line), and gas density (dashed line) 
along the axis y/H = 0. (b): Two-dimensional map of the index ip (see §3.3|) . 
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Fig. 13. — Height-time relation of the top of the flux tube. Dashed, solid, and dash-dotted 
lines represent cases for -B tu be = 67£> , 50-B , and 33B , respectively. Dotted line indicates 
the photospheric level. 
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Table 1. Summary of Cases 



Case 


B tuhc (300 G) a 


i? tubc (200 km) b 


q (0.005 km- 1 ) c 


1 


50 


5 


0.10 


2 


50 


5 


0.20 


3 


50 


5 


0.05 


4 


50 


5 


0.01 


5 


67 


5 


0.10 


6 


33 


5 


0.10 



a Magnetic field strength at the tube center. 
b Tube radius. 
c Twist parameter. 
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Table 2. Comparison of the Two Types of the Flux Emergence 



Undular Emergence Non-equilibrium Emergence 



Characteristic Variables 


of a Flux Sheet a 


of a Twisted Flux Tube b 


Initial Field Strength 


1.0 x 10 4 G 


1.5 x 10 4 G 


Total Magnetic Flux 


1.0 x 10 21 Mx c 


4.7 x 10 20 Mx 


Sheet Thickness / Tube Radius 


1000 km 


1000 km 


Initial Depth 


-20, 000 km 


-20, 000 km 


Arrival Time at the Surface 


4.9 x 10 4 s 


2.0 x 10 4 s 


Deceleration Depth 


-10,000 km 


-5000 km 



a Typical case of Paper I. 

b Typical case of this paper (twist q = 0.1/ H ). 

The flux sheet is assumed as a prism with a base 1000 km x 10,000 km 



